import math
from pyproj import Transformer, CRS
import csv
from shapely.geometry import Polygon, LineString, MultiLineString



# -----------------------------
# 坐标转换：WGS84 <-> 本地平面坐标 (Azimuthal Equidistant)
# -----------------------------
class CoordTransformer:
    def __init__(self, origin_point):
        lon0, lat0 = origin_point
        if abs(lat0) > 90:
            lon0, lat0 = lat0, lon0
        crs_wgs84 = CRS.from_epsg(4326)
        proj_str = f"+proj=aeqd +lat_0={lat0} +lon_0={lon0} +units=m +datum=WGS84"
        crs_local = CRS.from_proj4(proj_str)

        self.to_local_tf = Transformer.from_crs(crs_wgs84, crs_local, always_xy=True)
        self.to_wgs_tf = Transformer.from_crs(crs_local, crs_wgs84, always_xy=True)

    def to_local(self, lon, lat):
        x, y = self.to_local_tf.transform(lon, lat)
        return [x, y]

    def to_wgs(self, x, y):
        lon, lat = self.to_wgs_tf.transform(x, y)
        return [lon, lat]


# -----------------------------
# 几何运算函数
# -----------------------------
def slope_ax_by_c(p1, p2):
    """计算两点确定的直线方程 ax + by + c = 0"""
    if abs(p2[0] - p1[0]) < 1e-9:  # 垂直线
        return [None, 1, 0, -p1[0]]
    if abs(p2[1] - p1[1]) < 1e-9:  # 水平线
        return [0, 0, 1, -p1[1]]
    slope = (p2[1] - p1[1]) / (p2[0] - p1[0])
    a = slope
    b = -1
    c = p1[1] - slope * p1[0]
    return [slope, a, b, c]


def line_intersection(l1, l2):
    """计算两条直线的交点"""
    a1, b1, c1 = l1[1], l1[2], l1[3]
    a2, b2, c2 = l2[1], l2[2], l2[3]

    # 处理垂直线和水平线的特殊情况
    if a1 == 0 and b1 == 1:  # l1是水平线 y = -c1
        if a2 == 0 and b2 == 1:  # 两条水平线
            return None
        x = (-c2 - b2 * (-c1)) / a2 if a2 != 0 else None
        return [x, -c1] if x is not None else None
    elif b1 == 0:  # l1是垂直线 x = -c1
        if b2 == 0:  # 两条垂直线
            return None
        y = (-c2 - a2 * (-c1)) / b2 if b2 != 0 else None
        return [-c1, y] if y is not None else None

    # 一般情况
    D = a1 * b2 - a2 * b1
    if abs(D) < 1e-9:
        return None
    x = (b1 * c2 - b2 * c1) / D
    y = (a2 * c1 - a1 * c2) / D
    return [x, y]


def point_to_line_dist(a, b, c, x, y):
    """计算点到直线的距离"""
    return abs(a * x + b * y + c) / math.hypot(a, b)


def is_point_on_segment(p1, p2, pt, tol=1e-6):
    """检查点是否在线段上"""
    x, y = pt
    x1, y1 = p1
    x2, y2 = p2

    # 检查点是否在由p1和p2定义的矩形内
    if (min(x1, x2) - tol <= x <= max(x1, x2) + tol and
            min(y1, y2) - tol <= y <= max(y1, y2) + tol):
        # 检查共线性
        cross = (x - x1) * (y2 - y1) - (y - y1) * (x2 - x1)
        if abs(cross) < tol:
            return True
    return False


def is_inside_polygon(poly, pt):
    """检查点是否在多边形内部"""
    x, y = pt
    cnt = 0
    n = len(poly)
    for i in range(n):
        x1, y1 = poly[i]
        x2, y2 = poly[(i + 1) % n]
        if ((y1 > y) != (y2 > y)):
            xi = x1 + (y - y1) * (x2 - x1) / (y2 - y1)
            if xi > x:
                cnt += 1
    return cnt % 2 == 1


def offset_polygon(polygon_points, offset_distance):
    """使用shapely对多边形进行缩放"""
    polygon = Polygon(polygon_points)
    if not polygon.is_valid:
        polygon = polygon.buffer(0)

    # 计算缩放后的多边形
    if offset_distance > 0:
        offset_poly = polygon.buffer(offset_distance, join_style=2)
    elif offset_distance < 0:
        offset_poly = polygon.buffer(offset_distance, join_style=2)
    else:
        return polygon_points

    # 获取缩放后的多边形顶点
    if offset_poly.geom_type == 'Polygon':
        return list(offset_poly.exterior.coords)
    elif offset_poly.geom_type == 'MultiPolygon':
        # 取最大的多边形
        largest_poly = max(offset_poly.geoms, key=lambda p: p.area)
        return list(largest_poly.exterior.coords)
    else:
        return polygon_points


# -----------------------------
# 主流程：生成覆盖路径并保存到 CSV
# -----------------------------
def generate_coverage_path(params,wgs84_points):
    """
    生成覆盖路径

    参数:
    wgs84_points -- 多边形顶点坐标列表 [[lon, lat], ...]
    step_size -- 扫描线间距（米）
    offset_distance -- 多边形缩放距离（正值为外扩，负值为内缩）
    output_csv -- 输出CSV文件路径
    """
    step_size=params.get("step_size",3.0)
    offset_distance=params.get("offset_distance",0.0)
    if step_size == 0.0:
        step_size =3.0
    # 1. 坐标转换
    transformer = CoordTransformer(wgs84_points[0])
    local_pts = [transformer.to_local(lon, lat) for lon, lat in wgs84_points]

    # 2. 多边形缩放
    if offset_distance != 0:
        local_pts = offset_polygon(local_pts, offset_distance)
        # 移除最后一个重复的点（shapely返回的闭合多边形）
        if len(local_pts) > 1 and local_pts[0] == local_pts[-1]:
            local_pts = local_pts[:-1]

    n = len(local_pts)
    if n < 3:
        raise ValueError("缩放后的多边形顶点数不足，请减小offset_distance")

    # 3. 扫描方向：最长边
    max_d, idx1 = 0, 0
    for i in range(n):
        d = math.hypot(local_pts[(i + 1) % n][0] - local_pts[i][0],
                       local_pts[(i + 1) % n][1] - local_pts[i][1])
        if d > max_d:
            max_d, idx1 = d, i
    idx2 = (idx1 + 1) % n
    slope0, a, b, c0 = slope_ax_by_c(local_pts[idx1], local_pts[idx2])

    # 4. 终止截距
    others = [i for i in range(n) if i not in (idx1, idx2)]
    far_idx = max(others, key=lambda k: point_to_line_dist(a, b, c0, *local_pts[k]))
    xf, yf = local_pts[far_idx]
    stopping_c = -(a * xf + b * yf)

    # 5. 生成平行扫描线，调整步长
    lines = []
    delta_c = step_size * math.hypot(a, b)
    sign = 1 if stopping_c > c0 else -1
    c = c0 + sign * (delta_c * 0.5)
    while (sign == 1 and c < stopping_c) or (sign == -1 and c > stopping_c):
        lines.append([slope0, a, b, c])
        c += sign * delta_c

    # 6. 准备路径累积
    path_local = []
    dir_vec = (b, -a)  # 扫描方向向量

    # 7. 改进的路径生成算法
    for idx, scan in enumerate(lines):
        seg_pts = []

        # 计算扫描线与多边形边的交点
        for i in range(n):
            p1 = local_pts[i]
            p2 = local_pts[(i + 1) % n]
            edge = slope_ax_by_c(p1, p2)
            pt = line_intersection(scan, edge)

            if pt is not None:
                # 检查交点是否在边上
                if is_point_on_segment(p1, p2, pt):
                    # 检查是否已经包含该点（避免重复）
                    if not any(math.hypot(pt[0] - p[0], pt[1] - p[1]) < 1e-6 for p in seg_pts):
                        seg_pts.append(pt)

        # 如果没有交点，尝试更宽松的容差
        if len(seg_pts) < 2:
            for i in range(n):
                p1 = local_pts[i]
                p2 = local_pts[(i + 1) % n]
                edge = slope_ax_by_c(p1, p2)
                pt = line_intersection(scan, edge)

                if pt is not None:
                    # 使用更宽松的容差检查
                    if is_point_on_segment(p1, p2, pt, tol=1e-4):
                        if not any(math.hypot(pt[0] - p[0], pt[1] - p[1]) < 1e-4 for p in seg_pts):
                            seg_pts.append(pt)

        # 如果还是没有足够交点，跳过这条扫描线
        if len(seg_pts) < 2:
            print(f"警告：扫描线 {idx + 1} 只有 {len(seg_pts)} 个交点，跳过")
            continue

        # 确保交点数量为偶数（可能需要添加顶点）
        if len(seg_pts) % 2 != 0:
            # 找出最接近扫描线的顶点
            closest_pt = min(local_pts, key=lambda p: abs(a * p[0] + b * p[1] + scan[3]))
            seg_pts.append(closest_pt)

        # 投影后排序
        projs = [pt[0] * dir_vec[0] + pt[1] * dir_vec[1] for pt in seg_pts]
        sorted_pts = [pt for _, pt in sorted(zip(projs, seg_pts), key=lambda x: x[0])]

        # 交替方向
        if idx % 2 == 0:
            path_local.extend(sorted_pts)
        else:
            path_local.extend(sorted_pts[::-1])

    # 8. 回 WGS84 坐标
    path_wgs84 = [transformer.to_wgs(x, y) for x, y in path_local]

    # 9. 保存路径到 CSV
    # with open(output_csv, 'w', newline='') as f:
    #     writer = csv.writer(f)
    #     writer.writerow(["Longitude", "Latitude"])
    #     writer.writerows(path_wgs84)

    return path_wgs84
